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We present a general procedure for obtaining progressively more accurate functional expressions 
for the electron self-energy by iterative solution of Hedin's coupled equations. The iterative process 
starting from Hartree theory, which gives rise to the GW approximation, is continued further, and 
an explicit formula for the vertex function from the second full cycle is given. Calculated excitation 
energies for a Hubbard Hamiltonian demonstrate the convergence of the iterative process and provide 
further strong justification for the GW approximation. 



Much of the modern theory of many-body effects in 
the electronic structure of sohds relies on a closed set of 
coupled integral equations known as Hedin's equations 
0, which connect the Green's function G of a system 
of interacting electrons with the self-energy operator E, 
the polarization propagator P, the dynamically screened 
Coulomb interaction W , and a vertex function F. Si- 
multaneously solving Hedin's equations for a specified 
external potential in principle yields the exact Green's 
function and quasiparticle excitation spectrum without 
the need of actually calculating the many-electron wave- 
function. Unfortunately, however, the relation between 
the quantities is not just purely numerical but involves 
nontrivial functional derivatives, so that an automated 
numerical solution is not feasible and approximate func- 
tional expressions with simpler dependences have to be 
considered instead. Most calculations of quasiparticle ex- 
citations in real materials employ the GW approximation 
[01 , which uses intermediate operators from the first cycle 
of an iterative solution of Hedin's equations starting from 
Hartree theory as a zeroth order approximation. The 
GW approximation neglects diagrammatic vertex correc- 
tions both in the polarization propagator and the self- 
energy. Its theoretical foundation lies in the assumption 
that sufficient convergence has been reached after the ini- 
tial cycle, but rigorous evidence has so far been prevented 
by the inherent mathematical difficulties associated with 
a continuation of the iterative process. Some corrobo- 
ration stems from the surprisingly good agreement with 
experimental spectra for a wide range of semiconductors 
[H,D and simple metals [Q, while more recent studies of 
transition metals and their oxides highlighted deficiencies 
in the calculated spectra that clearly indicate a lack of 
convergence for materials with strong electronic correla- 
tion 1^. The GW approximation has since often been 
reinterpreted as the first order term in an expansion of 
the exact self-energy in terms of the screened Coulomb 
interaction, and considerable effort has been spent to 
include second order contributions. However, such at- 
tempts have failed to produce a general improvement in 



numerical accuracy due to far-reaching cancellation be- 
tween the additional terms 1^. In this Letter, we return 
to the original spirit of solving Hedin's equations itera- 
tively and present a scheme that generates progressively 
accurate functional expressions at arbitrary levels of iter- 
ation. Starting from Hartree theory, we obtain an explicit 
formula for the vertex function from the second iterative 
cycle, which mixes certain diagrams of different orders in 
the screened interaction. As exemplified for a Hubbard 
Hamiltonian, this approach indeed yields convergent ex- 
citation energies beyond the GW approximation. 

In our notation, the initial zeroth order self-energy 
and corresponding Green's function are labelled and 
G*-"-*, respectively. The sequence in which the operators 
are obtained during the first iterative cycle then starts 
with r(i), followed by P^'^\ W^^\ S^^), and finally G^^). 
For two reasons the principal mathematical difficulty in 
continuing this process towards convergence lies in the 
calculation of the vertex function. First, it is defined 
only implicitly through the Bcthe-Salpeter equation 

xG(")(7, 5)F("+i)(6, 7; 3) d(4, 5, 6, 7) (1) 

with the labels 1,2,... each denoting a set of position, 
time, and spin variables. While integral equations for 
other operators are readily solved by matrix inversion 
in Fourier space, the convolutions in (P cannot easily 
be disentangled, so that the computational expense is 
prohibitive. Second, it contains the functional deriva- 
tive 5Y.'^^^ / 5G'^^\ which is nontrivial because the Green's 
function G'"-* is not explicitly contained in E'"^ but only 
calculated from it by means of Dyson's equation 

G(")(8,9) = G(")(8,9) + Jg("H8,1) 

xAE(")(l,2)G(")(2,9)d(l,2), (2) 

where AE^"^ = E'"^ — E^°). A numerical treatment of 
the vertex function becomes feasible only if the functional 
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FIG. 1. The diagrammatic series (a) for SG^^^^ /5G'"^ can 
be summed using an expansion of Dyson's equation in the 
form (b) and has the exphcit solution (c) . By adding the first 
two terms on the right-hand side of (d) and applying the chain 
rule to the third, we obtain an integral equation (e) that can 
be solved by relation to the Bethe-Salpeter equation for n — 0. 
Single and double line arrows indicate the Green's functions 
G'"' and G^"\ respectively, the encircled E represents the 
self-energy correction AS^"^ = E'"' — and the shaded 
semicircle the vertex function r^""*"^-*. Functional derivatives 
are labelled explicitly. 

derivative is evaluated analytically, and if (0) can be 
solved explicitly for r'^"+^^ In the following we present 
a transformation that satisfies both requirements. 

All operators relevant in this context contain as their 
basic building block the zeroth order Green's function 
qW _ We therefore start by applying the chain rule 



5E(")(1,2) _ /■(5S(")(1,2) ,5G("'(8,9) 



<5G(") (4, 5) J 5G(") (8, 9) (5G(«) (4, 5) 



d(8,9). (3) 



The first functional derivative can in principle be evalu- 
ated at any level of iteration, so we focus on the second 
term. We perform the derivative of (||) with respect to 
G'"', which yields an integral equation for the four-point 
operator SG^^^ / 6G^"'^ as shown in diagrammatic form in 
Fig. |l|(a) . This series can be summed using an expansion 
of Dyson's equation in terms of G*^"-* with alternating 
signs [Fig. |l](b)]. The explicit solution is given in Fig. 
|l|(c). Employing Dyson's equation once more, we can 
write this relation more concisely as 

j|!^G(")(4,6)G(")(7,5)d(4,5) 
G(°) (8, 6)G(°) (7, 9) - / G(") (8, 1)G(°) (2, 9) 



lGH(Hr^^"^ (4, 6)G(") (7, 5) 2, 4, 5), (4) 

which still contains functional derivatives of the self- 
energies. Next we resubstitute AS^") = S^") - S^") and 



perform the integral with the vertex function r("+-'^) as 
required by the Bethe-Salpeter equation. These two steps 
yield the equation shown diagrammatically in Fig. |l|(d). 
We simplify the sum of the first two terms on the right- 
hand side using the relation (|^) and rewrite the functional 
derivative in the third by applying the chain rule 

^SW(1,2) _ /■ JS](")(1,2) ^GW(8,9) 

,5G(")(4,5) 7,5G(o)(8,9)5G(")(4,5) ^ ' ^ ' 

In this way we finally obtain the integral equation shown 
in Fig. |l|(e), which is closely related to the Bethe-Salpeter 
equation for n — 0. By comparison, we find 

G(") (8, 6)G(°) (7, 9)r(i) (6, 7; 3) d(6, 7) 



(5GW(8,9) („) Q)Qin)(j 
JG(")(4,5) ^ ' ^ ^ ' ^ 
X d(4,5,6,7). 



n+l] 



(6, 7; 3) 



(6) 



We can now use this identity to sum the diagrammatic 
series in (|l|) and rewrite it in the alternative form 

r(«+i) (1, 2; 3) = 6{1,2)S{1, 3) + J^^^^^G^o) g) 

xG(°) (7, 5)r(i)(6, 7; 3) d(4, 5, 6, 7). (7) 

This expression is remarkably similar to the original in- 
tegral equation, except that all operators but S*^"^ on the 
right-hand side are replaced by their lowest order equiv- 
alents. While successive iterations dress each propagator 
with new sets of diagrams, the identity (|^) indicates far- 
reaching cancellation between the expansion terms from 
individual propagators, yielding a much simpler expres- 
sion for the vertex corrections than originally anticipated. 
The emergence of propagators from the lowest iterative 
cycle reduces the numerical expense substantially, as in 
practice mean-field theories are used as a zeroth order ap- 
proximation. In this case G^°^ contains no satellite spec- 
trum but only a set of robust quasiparticle excitations. 
The transformation also satisfies our requirements for a 
numerical treatment by giving an explicit definition for 
the vertex function and expressing the functional deriva- 
tive in a way that can be evaluated at higher orders. In 
the self-consistency limit n — > 00, (Q) implies a relation 
between the exact self-energy and vertex function. 

In condensed matter physics, iteration conventionally 
starts from Hartree theory as a zeroth order approxi- 
mation, so that = and G(°) = G". Due to the 
vanishing self-energy the functional derivative in (|l|) is 
identically zero, yielding a trivial vertex function. The 
subsequent iteration generates the GW approximation 

P(i)(l,2) = -iG(")(1,2)G(°)(2,1) 



H^(i)(l,2) = i;(l,2)+ /m^(i)(1,3) 



xp(i)(3,4)v(4,2)rf(3,4) 
I](i)(l,2) = iG(")(l,2)iy(i)(l+,2), 



(8) 
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FIG. 2. Vertex correction F'^-* to the GW self-energy from 
the second iterative cycle of Hedin's equations. Single line 
arrows and wiggly lines represent the zeroth order Green's 
function G^'^^ and the screened Coulomb interaction W^^^ in 
the random phase approximation, respectively. 



where in the last equation 1+ implies that a positive in- 
finitesimal is added to the time variable, v is the bare 
Coulomb interaction. At this level the self-energy is mod- 
elled as the product of the Green's function G*^"^ and 
the screened interaction W^^"^ in the random phase ap- 
proximation (RPA). The functional form is reminiscent 
of the Fock potential, but the electronic exchange in- 
cludes dynamic screening and so reaches beyond the lim- 
its of mean-field theories. Using this definition of the 
self-energy, it is easy to verify that its derivative with 
respect to the zeroth order Green's function is given by 

~;.,.(MM(2,5„W..(1,2) 

+ T4^(i) (1,4)14^(1^5, 2)1 G(°)(5,4). (11) 

The corresponding vertex function F^^^ that underlies the 
second iterative cycle is shown in diagrammatic form in 
Fig. H This finite set of vertex corrections is distinct 
from that obtained through expansion by orders of the 
screened interaction in that it comprises selected terms 
which are of zeroth, first, and second order in FF^^^. 

Although the self-energy ( |lO| ) was obtained by itera- 
tion starting from Hartree theory, in practice it is more 
often evaluated using a zeroth order Green's function 
Qio) — G^^""" from a previous density-functional calcula- 
tion with S'^*'^ = V-^'~" equal to the exchange-correlation 
potential. Although this substitution violates the orig- 
inal spirit of the iterative scheme, physical arguments 
suggest only a small deviation to the propagators prop- 
erly derived from density-functional theory as a zeroth 
order approximation [^. Since the form of S*^^^ remains 
identical to that considered before, its derivative is still 
given by (pd]). The nontrivial first order vertex F^^^ re- 
quired to evaluate (0) was derived in Ref. and indeed 
found to have insignificant numerical impact on the band 
gap of silicon. It is probably also negligible in the sec- 
ond iteration of such a calculation, so that the expression 
shown in Fig. ^ may still be used for F^^^ . 

In the following we will supplement our discussion with 
the numerical investigation of a system of strongly corre- 
lated electrons that explores the effects of the vertex cor- 
rection F*^^). As the evaluation of higher order diagrams 
for real materials would demand computational resources 



beyond the scope of the present study, we consider a two- 
dimensional square array of 3 x 3 atoms described by the 
Hubbard Hamiltonian 



(R,R'),<T 



R 



(12) 



The model captures the essential physical features of 
materials with strong electronic correlation such as the 
transition metals, which are known to be inadequately 
described in the GW approximation but is simple 
enough to allow the calculation of full optical spectra 
rather than just individual matrix elements as in Refs. 

Here, and cro- are the creation and annihi- 
lation operators for an electron at site R with spin cr, 
and we define fiji^ — c^^cro-. The notation (R, R') im- 
plies a sum over nearest neighbors only. This model was 
introduced in Ref. j|] to compare variations of the GW 
approximation, and its performance within the first iter- 
ative cycle of Hedin's equations is thus well understood. 
The single band of the cluster can accomodate up to 18 
electrons; we consider a system of 16 electrons. The high 
fractional band filling resembles that of the d orbitals in 
the late transition metals. Although we use open bound- 
ary conditions, the on-site energies er are chosen in such 
a way as to yield uniform occupation numbers in the 
Hartree approximation, as expected in infinite systems. 
For reference, the on-site energy is 2t for corner sites, t 
on edge sites, and in the center of the cluster. We use 
medium interaction U = At and set t = 1. 

In Fig. ^ we compare the exact screened interaction 
W with its RPA counterpart from the first iteration of 
Hedin's equations, in which the polarization propagator 
takes the form (^), and with the second iteration. Here 

P(')(I,2) = -I /G(i)(2,3)G(i)(4,2)F(2)(3,4;I)d(3,4). 



(13) 

As in Ref. ||], G^^^ was obtained from a shifted zeroth 
order Green's function in order to avoid problems with 
differing chemical potentials in the perturbation series. 
The figure shows the diagonal matrix element for the 
central site of the cluster, which because of the chosen 
geometry corresponds most closely to the screening in 
extended systems. Other matrix elements exhibit similar 
behavior, and the displayed curves are thus representa- 
tive. For comparison, we also show results from an ex- 
pansion of the polarization propagator beyond the RPA 
to first order in VF^^^. This approach is qualitatively dis- 
tinct in that only the vertex fmiction becomes dressed, 
while in the iterative scheme both F and the internal 
propagators G in ( p^ are simultaneously updated. 

The exact screening is dominated by a pair of strong 
plasmon peaks at 3 and just under 6 eV, but three satel- 
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FIG. 3. The exact screened interaction compared to the 
RPA, an expansion of the polarization to first order in W, and 
the second iteration of Hedin's equations. The latter yields 
more accurate plasmon energies and a qualitatively correct 
satellite spectrum. Inset: The shift in plasmon energies is due 
to the explicit vertex function and remains when the diagrams 
are evaluated with G'-"-' rather than G^^\ while the satellite 
peaks result from internal dressed propagators. 

lites at higher energy can be identified. While qualita- 
tively acceptable, the RPA as a first approximation ig- 
nores the satellite spectrum and places the plasmons too 
high by about 1 eV. The latter deficiency is somewhat 
improved by the further expansion of the polarization 
propagator, but the description of the satellites remains 
poor and is not even qualitatively correct. This is in line 
with previous observations of far-reaching cancellation 
between the additional terms ||^. 

In comparison, the second iteration is more effective in 
shifting the plasmon energies, particularly for the lower 
peak, and also yields a better low-frequency limit for the 
static dielectric function. Furthermore, we observe the 
emergence of a satellite spectrum that is in good agree- 
ment with the exact curve concerning the number and 
position of features. Their exaggerated spectral weight 
can be traced to the satellites in the GW Green's func- 
tion, which for this model are overestimated by a similar 
factor ||]. To isolate the effect of the vertex function, 
we have also evaluated (Q3h using G^o) rather than G'^^\ 
The resultant dielectric function, displayed in the inset 
in Fig. ^ retains the full shift in plasmon energies but 
no satellites, which are due to the internal dressed prop- 
agators. Incidentally, the graph also shows part of the 
apparent plasmon spectral weight to have moved across 
the real axis. This incorrect analytic structure is an un- 
fortunate but common consequence of nontrivial vertex 
corrections and due to the occurrence of higher order 
poles in P^'^l The same effect is well documented for the 
expansion in terms of the screened Coulomb interaction 



where it is also observed in the present calculation. 
However, we have confirmed that the integral of the spec- 
tral function — ImW^'^^^(tj)/7r over the positive half-axis 
is correctly greater than zero for all diagonal matrix el- 
ements, and as intervals of negative spectral weight are 
few and embedded between features with correct analytic 
behavior, they do not dominate convolutions in the later 
course of the iteration. 

In summary, we have presented a scheme for the sys- 
tematic construction of vertex corrections by iterative so- 
lution of Hedin's coupled equations, and we have given 
explicit formulae for the propagators beyond the GW ap- 
proximation. Numerical results for a model of strongly 
correlated electrons indicate that this method not only 
yields improved excitation energies but is also more pow- 
erful than a comparable expansion by orders of the 
screened Coulomb interaction, in particular by gener- 
ating a superior satellite spectrum. On a fundamental 
level, these findings provide the first direct evidence of 
convergence of the iterative approach and thus give fur- 
ther theoretical justification for the GW approximation. 
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